Multiscale Modeling of a Nanoelectromechanical Shuttle 
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In this article, we report a theoretical analysis of a nanoelectromechanical shuttle based on a 
multiscale model that combines microscopic electronic structure data with macroscopic dynamics. 
The microscopic part utilizes a (static) density functional description to obtain the energy levels and 
orbitals of the shuttling particle together with the forces acting on the particle. The macroscopic 
part combines stochastic charge dynamics that incorporates the microscopically evaluated tunneling 
rates with a Newtonian dynamics. 

We have applied the multiscale model to describe the shuttling of a single copper atom between 
two gold-like jellium electrodes. We find that energy spectrum and particle surface interaction 
greatly influence shuttling dynamics; in the specific example that we studied the shuttling is found 
to involve only charge states Q = and Q — +e. The system is found to exhibit two quasi-stable 
shuttling modes, a fundamental one and an excited one with a larger amplitude of mechanical 
motion, with random transitions between them. 



I. INTRODUCTION 

Nanoelectromechanical systems that combine electri- 
cal and mechanical functionalities on the nanometer scale 
have in the recent years attracted a great deal of theo- 
retical and experimental interest The nanoelectrome- 
chanical shuttle is a structure that resembles a single elec- 
tron transistor but incorporates mechanical motion of the 
central island. Previous theoretical works on the shuttle 
have shown that in the presence of an DC applied bias the 
charge and velocity of the central island are correlated, 
Q(t)Z(t) ^ 0, which implies that the shuttle absorbs en- 
ergy from the DC field and converts it into mechanical 
motion. The shuttle motion facilitates charge transfer 
through the system, and signatures of mechanical mo- 
tion can be seen both in the current-voltage characteris- 
tics and in the noise properties of the device^*a 

Several theoretical studies have been carried out for 
different setups of the shuttle since the first description 
of this phenomenon. The theoretical studies cover dif- 
ferent size regimes of the shuttle, featuring coherent or 
sequential^ tunneling and quantum mechanically«2iifl or 
classicallyii*ia described mechanical motion. The studies 
have shown that the shuttle instability strongly depends 
on the bias voltage and the system setup. This sensitiv- 
ity also renders the shuttle behavior dependent on the 
precise description of the problem. 

Experimental evidence of coupling between vibrational 
degrees of freedom and electron transfer has been found 
for both microscopioiiiii and macroscopic*^ systems. In 
particular, the experiment by Park et at, reference Il6t 
using a Cgo molecule between gold electrodes has demon- 
strated the type of coupling that has been considered by 
many theoretical studies and has increased the interest 
for a molecular shuttle pi*2*i2*i2*i£ 

In the shuttle geometry the mechanical motion is on a 
nearly macroscopic time scale, typically from picoseconds 
for small molecules to nanoseconds for large molecules 
such as carbon nanotubes. The motion is excited due to 
tunneling events between the mobile object and the sta- 



tionary electrodes, which have a typical timescale of fem- 
toseconds and are determined by the electronic structures 
of the mobile molecule and the electrodes. Hence, a the- 
oretical description of the shuttle system naturally calls 
for multiscale methods that combine the fast electronic 
time scales with the slower mechanical ones. Thus far 
research has concentrated on the slow degrees of freedom 
while dealing with the fast ones in a phenomenological 
approximation. 

The two main issues addressed in this work are the im- 
pact of the electronic structure of the central island on 
the shuttling motion, and an analysis of the short range 
interactions between the central island and the station- 
ary electrodes. The first issue we will address by describ- 
ing the central island using density functional theory^iS 
which provides information on the energy spectrum of 
the island as well as structure of the relevant orbitals. 
The interaction between the island and the electrodes is 
described in part by a phenomenological Born-Mayer po- 
tential combined with image charge effects, and in part by 
a model that transfers mechanical energy from the shuttle 
to lattice vibrations in the electrodes, thereby dissipat- 
ing energy of the shuttle system and preventing catas- 
trophic runaway. Due to the phenomenological descrip- 
tion of surface interactions, some physical effects such as 
chemisorption are not properly accounted for which lim- 
its the applicability of the present model to materials for 
which chemisorption can be neglected. A way to over- 
come this problem in a future study would be to incor- 
porate a time-dependent DFT 20 module that describes 
both the island an the electrode during the crucial parts 
of the shuttle cycle; at present, however, that type of de- 
scription is prohibitively expensive from a computational 
point of view. 

The DFT data is used to evaluate tunneling matrix 
elements and tunneling rates between the central island 
and the electrodes. These are then inserted to a dynamic 
module that describes both the charge dynamics in terms 
of stochastic tunneling events and mechanical motion us- 
ing molecular dynamics. The resulting macroscopic dy- 
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FIG. 1: The system. A copper atom is placed in between two 
jellium surfaces 15 A apart. A buffer with the width of 2.5 
A is added on both side of the gap for the DFT calculations 
in order to localize the electrons. A small smoothing area 
between gap and barrier is needed to speed up calculations. 



namics is implemented as a dynamic Monte Carlo algo- 
rithm that uses the output of a series of static DFT cal- 
culations as its input. This stochastic dynamics can be 
used to address both average transport properties such 
as current and random fluctuations, or noise, which both 
exhibit clear signatures of shuttling according to phe- 
nomenological theories. 



II. METHOD 



A. Setup 



For computational efficiency, we have chosen to fo- 
cus on the simplest possible system where the central 
island comprises just one atom. However, the methods 
and qualitative results should be applicable also to more 
complex system. The system we consider includes two 
electrodes 15 A apart, described as semi-infinite jellium 
slabs. The central island is a copper atom that can move 
in a direction normal to the electrode surfaces. For the 
electrodes, the Wigner-Seitz radius is set to 3 a.u. and 
the electrode work function W is set to 3.5 eV—i The 
fermi energy, e^, is calculated from the Wigner-Seitz ra- 
dius while other material specific electrode parameters 
are taken from gold. A bias voltage of 3 V is applied 
over the gap, the potential dropping from left to right. 
The region described by the DFT module consists of the 
region between the two electrodes plus a buffer region of 
2.5 A inside each of the electrodes as show in Fig. ^ 
The buffer regions are needed so that the plane-wave- 
based code can better describe orbitals that are localized 
in the inter-electrode gap. 



B. Electronic structure calculations 

The electrodes and the space between them are treated 
separately: the electronic structure of the shuttling ob- 
ject is described using a density functional code that is 
limited to the inter-electrode space while the electrodes 
themselves are described analytically within the jellium 
approximation. For the electronic structure calculation 
we use the DaCapo code 2 ^ with the PW91 exchange cor- 
relation functional 2 — and employ the adiabatic approx- 
imation to separate the electronic structure calculation 
from the dynamics description. 

The effect of the electrodes on the island is modeled 
as a one-electron potential and inserted directly into the 
effective potential in the Kohn-Sham equations. The po- 
tential is divided into two parts: the interaction between 
an electron and a metal surface, and the interaction be- 
tween an electron and the induced charge arising from the 
remaining charges in the space between the electrodes. 
For the former part we use the saturated image barrier 



Vj(z) 



167reo(z 
V 



zo) 



-\(z-zO)\ 



A e B{z-z ) + 1 



Z > Zq 



Z < Zq, 



(1) 

suggested by Jones and Jennings H Here the param- 
eters are related to the work function W and Fermi 
energy tp of the electrodes through Vq = W + ep, 
A = -1 + 167re Vo/<7 2 A and B = 8ne Vo/q 2 A where the 
parameters zq and A are found by fitting eq. with 
DFT calculated tabulated data on the effective poten- 
tial, v e gm^ The values were determined to A = 1.6 A -1 
and z = 0.34 A outside the jellium background Zf,. The 
suppression factor (1 — exp(— A (z — zq))), which multi- 
plies the classical image potential, tends towards zero as 
the mirrored charge closes in on the surface. The image 
plane, zq, is used as the effective surface of the material, 
in compliance with Lang and Kohn— — The value of zq 
derived for eq. Q differs from the value derived by Lang 
and Kohn for e.g. image charge potential. 27 However, for 
the interactions with the induced charge, the individual 
contributions from the electrons and the nucleus should 
cancel out far from the surface and it is therefore suitable 
to keep only one parameter for the effective surface. 

For the second term, the potential at r due to a charge 
at position r', we use a form i>(r;r') that satisfies Pois- 
son's equation in the region outside the electrodes and 
reduces to the saturated image potential if the mirrored 
charge is at the same point where the potential is mea- 
sured giving 
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The resulting potential for interaction with the core is 
V In (r)=Q v v(r;Z), (3) 
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where Q v is the core charge of the island, i.e., the number 
of protons minus the number of core electrons. 

The model for the interaction with the induced charge 
from the other electrons is obtained similarly. Here we 
assume that the electron distribution of the island is well 
approximated with a Gaussian, as is typically the case, 
and integrate over the mirrored charge to obtain 

V Ie (r) = -(Q e - 1) J dr'v(r; r')p(r'). (4) 

where 

*'> = w- 4,z - w < 5 » 

is a form function and Q e is the number of valence elec- 
trons for the central island. The width a(Q) is calculated 
by fitting a Gaussian to the unperturbed electron den- 
sity for different charge states, Q denoting the number 
of extra charge units on the central island. Assuming 
localization of electrons to the island or the leads (no 
chemisorption) , Q is strictly integer. 
The added bias voltage is 

V Has (z) = E(z - A gap /2) (6) 

where the electric field is E = V/(A gap ) and A gap the dis- 
tance between the electrodes. Finally, the Pauli repulsion 
that confines electrons within the gap is for the electronic 
structure calculations described as a repulsive square po- 
tential wall placed at Zf,. The repulsion term is important 
for the separation of electronic structure calculation and 
charge transfer mediated by tunneling, however, the form 
of the repulsion is less important as most of the tunneling 
events take place when the island is relatively far from 
the surface on an atomic scale. 

Owing to the time-consuming and time-independent 
character of DFT, it is not possible to determine the sys- 
tem properties continuously. Instead, the simulations are 
performed on a number of positions and charge states. A 
continuous description is produced with interpolation. 



C. Forces 

The forces on the atomic core are calculated using the 
DFT code. In the Kohn-Sham single-particle equations 
all charges but one electron are treated as a mean-field 
static charge distribution. We implement the potential 
due to surface interactions as an external field. This re- 
sults in proper (mean field) orbitals and energy eigenval- 
ues but incorrect forces and total energies. Therefore, 
a correction term in the form of dV g£ T — is added 
to the forces calculated by the electronic structure code. 
Here, Vbft is the induced charge potential used in the 
DFT calculation, the sum of eq. {Q, © and Q}. The ac- 
tual potential, V e i is obtained as a variational derivative 



of the full electronic energy as 

V el ~ J dV KZ; r') + v(r>; Z)] p(r') 
-Q v v{Z;Z)+g(r) 

where the last term g(r) is independent of Z and does 
not contribute to the forces. 

A smooth many-body short-range repulsion is added 
to the total force. We have chosen a Born-Mayer type 
pair potential 2 ^ and integrated over the electrode surface 
which yields the force 

F{z) = f Q e^+^cj/p) f^z/p + 2Trzpn s eS-V z2 + a2 /( 2 *yp^ . 

(8) 

For a Cu-atom outside a gold surface, the effective radii 
are ocu — 0.77 A and ctau — 1-37 A, n s — 2/ a 2 is the 
surface density of atoms and a; at = 4.078 A is the lat- 
tice constant for a fee [100] gold surface. 2 ?'?? The closest 
electrode lattice site is considered separately (in order to 
maximize the localization) while the other lattice sites 
are handled as a continuum. The values fo and p are 
6.95 10~ n N and 0.3 A respectively.? 1 !? 2 !??!? 4 



D. Transition rates 

Charge transfer rates between the central island and 
electrodes are calculated with the transfer Hamiltonian 
method^*^ using the overlap between the atomic and 
electrode wave functions. The transition rates to and 
from a specific atomic orbital are 

(9) 

where 

M pq = { d 3 r^ a (V a -V m )i> q , m - (10) 

J z£R gap 

Here, ipp,a is a Kohn-Sham orbital (localized on the is- 
land) with eigenenergy E p , ipq,m is a metal wave func- 
tion and n/ is the relevant number of states available for 
tunneling in the (often degenerate) orbital ip p , a - The po- 
tentials V a and V m are the effective potentials for ip p ^ a 
and ipq t m respectively. The atomic orbital ipp,a and the 
effective gap potential V a are extracted from the DFT 
simulations. The transition rates are calculated numeri- 
cally for all energetically allowed transitions determined 
by comparing the electrode chemical potentials with the 
Kohn-Sham eigenvalues. 

The electrode wave functions are calculated analyti- 
cally from a square potential with a finite barrier at the 
electrode surface and a hard wall at — oo. The electrode 
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wave functions in the z-direction are 



with k z = 
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offset between the physical and the geometrical surface 
is implemented to attain charge neutrality of the jellium 
slabs^i Periodic boundary conditions are assumed par- 
allel to the surface. 



E. Dynamics 

A dynamic Monte Carlo approach is used to calculate 
the shuttle dynamics. Input parameters to this module 
consist of core forces and transition rates as functions of 
the island position. The motion of the central island is 
described classically by 



mx = F ext (x, Q) + F dissip (x, Q) 



(12) 



where F ext are the core forces given by the previous cal- 
culations and F^^p is a dissipation term. The island 
position x(t) is calculated by numerically integrating the 
equation of motion, while the island charge Q{t) is al- 
lowed to change stochastically using the tunnel rates de- 
termined above. This results in a coupled stochastic dy- 
namics for the mechanical and electrical degrees of free- 
dom. 

As the shuttle absorbs energy from the bias voltage, 
the dissipation term is essential for the stability of the is- 
land motional Earlier theoretical work has mainly used 
viscous damping, —fix^i^^^J^ In this work a different 
model based on mechanical damping is used, coupling 
a simple model of the surface to the island equation of 
motion via the surface forces as 



mx = F tot (x) - X F> ot (x)G(X F{ ot (x)x) 
MX = -MX-kX-F tot (x-X), 



(13) 



where x and m are the position and mass of the central 
island, X and M are the position and mass of the surface, 
and k is an effective surface spring constant. The 0- 
function in the first equation restricts the energy flow 
so that the shuttle energy can be transferred to lattice 
vibrations of the electrodes (phonon emission) while the 
opposite process of phonon absorption is forbidden. 

The equations have been formulated in terms of the 
surface element nearest to the shuttling object, which 
implies that the effective mass of the surface depends 
on the separation between the mobile island and the sur- 
face; this can be determined by assuming an elastic model 
for the surface and requiring that (i) the instantaneous 
displacement X(t) agrees with that of the surface atom 
nearest to the shuttling object, and (ii) the total momen- 
tum of the surface is MX. The elastic parameters for 



surface atoms have been chosen to correspond to those 
of gold, and effective parameters in the simplified model 
are consequently M(X) = rriAuF{X) / F max (X), where 
F(X) is the total force between the shuttling object and 
the surface and F max (X) is the force on the surface atom 
nearest to the shuttle and k = k a u F {X) / F max {X) where 
kAu = SniAuVs /° 2 i s t ne spring constant obtained from 
the sound velocity v s and lattice constant a. 

The main qualitative difference between our dissipa- 
tion model and viscous damping is that in our model 
dissipation occurs primarily when the island-electrode in- 
teraction is strongest. This influences the threshold volt- 
age for onset of shuttling, and also renders the threshold 
dependent on the initial conditions. 



III. RESULTS 

The different parts of the one-electron potential are de- 
picted in Fig. J5J ■ The small widths of the form function 
((7 ~ 0.24-0.27 A for an unperturbed pseudopotential, 
Q = 1 ... — 1, and ~ 0.30-0.35 A for the double junction 
potential) imply that sufficiently far from the surface the 
point charge approximation would be quite accurate: for 
distances > 2 A from the surface, the spatial distribu- 
tion of charge has little effect on the potential. For is- 
land positions close to the electrode surfaces, the spread 
in the valence electron distribution and the rapid satura- 
tion effectively bares the core image making the effective 
potential strongly repellent. 

The resulting Kohn-Sham eigenvalues, depicted in Fig. 
EH are used as energy spectra for the central island. Com- 
parison between the eigenvalues an the electrode chem- 
ical potentials gives the possible transitions. Full relax- 
ation into the A^/2 lowest bands is assumed instanta- 
neous, where N is the number of valence electrons (11 
for the used Cu GGA pseudopotential, Q = 0). Higher 
bands are treated as excitations4i The temperature is 
taken to be zero in the treatment of tunneling events; 
however, in the DFT calculation a finite temperature is 
needed for convergence. 

A small correction is needed for some calculations in 
order to use the equilibrium DFT calculations within a 
dynamic picture. For some island positions and Q < 0, 
it is energetically favorable to place some of the extra 
charge in the surface potential well outside the positive 
electrode surface instead of on the central island. How- 
ever, the time scale for this direct equilibration between 
leads is very long. In order to find the energy spectra 
and orbitals that are relevant for the dynamic evolution, 
the surface well near the left electrode surface is manually 
suppressed for core positions near the right electrode; due 
to the polarity of the applied bias, similar problems do 
not arise for core positions near the left electrode. The 
possibility of transitions directly between the electrodes 
is kept, but the transition rates are small enough to be 
of no importance for the results. 

For the tunneling rate calculation the Kohn-Sham 
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FIG. 2: The Kohn-Sham one-electron effective potential com- 
prises several parts rendering diverse behaviors for different 
positions and charge states. The above figures are for Q — 
and Z = 10.17 A. The z-axis is the direction of island motion, 
r is parallel to the electrodes. In (a) and (b) z encompasses 
the gap while (c) includes the entire DFT cell with the 2.5 A 
buffer regions on both sides of the gap. (a) Equation Q and 
the bias voltage, (b) One electron interaction with induced 
charge due to other system charges, (c) The effective poten- 
tial as used by DaCapo. Close to the electrodes eq. Q forms 
deep wells. 
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FIG. 3: The lowest eight non-spin polarized Kohn-Sham 
eigenvalues for nine positions of the central island. Levels in- 
dicated by filled circles are fully occupied, empty circles half- 
occupied and crosses unoccupied. The upper and lower solid 
lines are fiR and [il respectively, (a) Q = 0. The lower lying 
excitations (core positions in the left half of the gap) have 
been calculated including the surface potential well near the 
electrode. The higher lying excitations (core positions near 
right electrode) have been calculated without the surface po- 
tential well near the left electrode as transitions directly from 
the right electrode to the left surface well have very small 
tunneling rates. For Z = 7.5 A and Z = 10.17 A results are 
depicted for both approaches (solid and dashed lines). For 
the positions closest to the electrodes, the width of the sur- 
face potential wells causes a substantial drop in eigenenergies. 
(b) Q = 1. 



eigenvalues are regarded as electron energies, which is 
known to be rigorously correct for the ionization poten- 
tial involving the HOMO level^i and believed to be rea- 
sonably accurate for the other levels as well4£ We assume 
that tunneling rates are sufficiently low so that the island 
fully relaxes between each tunneling event to the config- 
uration determined by time-independent DFT. The time 
scale for this relaxation is typically in the femtosecond 
range which is fast compared to the tunneling rates ex- 
cept for core positions very close to the electrode sur- 
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faces; however, since energetics severely limits the pos- 
sible tunneling processes, the instantaneous relaxation 
approximation is reasonably well justified for all core po- 
sitions. For the chosen bias voltage, the possible charge 
transitions for the island are 1 — *■ and — * 1. 

It is interesting to notice the asymmetry of attainable 
charge states that arises from the asymmetry of the en- 
ergy spectrum of the Cu atom: the dynamical evolution 
only involves charge states Q = and Q = 1 but not 
Q = — 1. The symmetric expression for charging energy 
E = Q 2 /2C used for larger metallic grains is only jus- 
tified if the level spacing on the island is small enough 
so that the electrostatic energy scales dominate. This 
asymmetry implies that the shuttling is asymmetric also 
in the sense that energy is absorbed from the DC field 
only during half a period which makes the system more 
sensitive to dissipative mechanisms. 

The occupied Kohn-Sham eigenfunctions are identi- 
fied as d- and s-orbitals in accordance with the expected 
electron configuration of Cu, 3d 10 4s 1 . Close to the elec- 
trode surfaces, the orbitals deform against the repulsion 
wall. For all but the closest position to the electrodes the 
4s-orbital gives the widest electron distribution and the 
largest contribution to the transition rates. 

The core forces for the central positions are strongly 
dependent on the derealization of the electron distri- 
bution of the central island (Fig. QJ. For the positive 
ion with Q = 1 the dominant force is the electrical bias 
while the potential of Q = —1 is a nearly symmetric im- 
age charge potential. For Q = 0, the sign depends on 
the description of the surface interactions, and with the 
interaction model we have chosen the neutral atom feels 
a slight net force towards the negatively charged right 
electrode. The repulsion from the surface is dominant for 
the two outermost positions on either side giving a ph- 
ysisorption minimum between 3-5 A from the electrode 
surface. 

The transition rates are much less sensitive to the sur- 
face description than the forces (Fig. [S]). Their distance 
dependence is approximately exponential as assumed by 
effective theories 43 with slight saturation for core posi- 
tions nearest to the electrodes with a tunneling length 
that is approximately 0.4 A with some variation for the 
different allowed transitions. Near the electrodes the en- 
ergetics considerations inhibit tunneling, as seen in Fig. 
[31 which can be viewed as a molecular equivalent of 
Coulomb blockade. 

In the dynamics simulations, the calculated forces and 
transition rates are joined. The result is indeed a stable 
shuttling regime where Q(t)Z(t) ^ 0. We have performed 
dynamical simulations starting from a variety of initial 
states, and seen that for most starting conditions the re- 
sults are quite similar: as a rule, the model does indeed 
shuttle electrons (Fig. [fJJ). However, for some initial con- 
figurations such as Z(t = 0) ~ 10 A and Q(t = 0) = the 
applied bias of 3 V is not sufficient to initiate shuttling. 
One of the more prominent differences between our re- 
sults and previous works is the complexity of the forces, 
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FIG. 4: Total forces on the grain core, center positions. Po- 
sitions not included Z = 0.75 A and Z = 14.25 A reach 
±(500 — 650) nN. Lines between positions correspond to the 
used interpolation scheme. 
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FIG. 5: Transition rates: electron current from right to 
left, (a) Squares: Transitions from negative right lead to is- 
land. Close to the negative lead, current mediating transfer is 
blocked by energetics (Coulomb blockade). Instead, an elec- 
tron can transfer against the bias back to the negative lead, 
(circles), (b) Electron transitions from island to positive left 
lead. 
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FIG. 6: Shuttle regime for initial conditions Z(0) — 0, Q(0) = 
and v(0) = 0. (a) Island position, (b) Charge state as a 
function of time. The shuttle carries one electron per period, 
(c) Island velocity. The main acceleration and deceleration 
are close to the electrodes. 



particularly near the electrode surfaces. Since the sur- 
face description we use is adapted from static analyzes, 
it is unclear how accurately it captures the interactions 
in a dynamic situation, and the detailed results are some- 
what uncertain. A slightly different potential renders the 
forces on the neutral atom positive over a larger range of 
positions, and the range of initial conditions that would 
result in shuttling would be smaller, implying that the 
threshold voltage for shuttling depends sensitively on the 
model for surface-island interactions and on the initial 
conditions. 

The detailed structure of the forces of the middle po- 
sitions is less important after shuttling is well estab- 
lished. The main forces become the close-range expo- 
nential forces of the electrodes and the applied electric 
field. For the asymmetric shuttle, energy is absorbed 
by the charged shuttle from the field during half a cycle 
while during the other half-cycle, after an clastic collision 
with the electron surface, a neutral shuttle moves nearly 
freely in the opposite direction. The energy loss during 
the shuttle-electrode collision cannot exceed the energy 
absorbed from the field if a stable periodic motion is to 
be established. 

The distribution of positions at which tunneling events 
take place depends on the transition rates, and the width 
of the distribution is connected to the spatial derivatives 
of the rates, i.e. on the tunneling lengths. The steeper 
1 — ► gives a more compact distribution as seen in Fig. 

13 

The shuttle reaches a stable shuttling motion quickly 
with a current of ~0.19 fiA and an amplitude of ~11.1 
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FIG. 7: The statistical distribution of charge transfer loca- 
tions. The solid and dashed line corresponds to — > 1 and 
1 — » transitions respectively. The quicker growing prob- 
ability of a 1 — > transition compresses the distribution of 
event locations. 



A, (see Fig. (SJ . The random character of the transition 
processes influences the turning points very little. There 
is, however, a possibility for the system to undergo semi- 
stable excitations due to randomness of transfer events 
and the position dependence of the energy spectrum near 
the right lead (Fig. |5J. Very close to the negative (right) 
lead, there is a possibility of a process in which an elec- 
tron first tunnels from the electrode to an initially pos- 
itively charged (Q = 1) shuttle that continues to move 
towards the right lead, followed by tunneling against the 
bias back into the lead, and finally a new tunneling event 
after the shuttle has changed its direction of motion. 
During the time that the charged shuttle spends near 
the electrode surface after the second tunneling event, it 
experiences a larger force than a neutral shuttle would, 
which allows it to absorb more energy from the poten- 
tial and results in an enhanced shuttling amplitude. The 
increase in amplitude enhances the possibility for this se- 
quence of three tunneling events to take place also in the 
next period. The excitation lasts until a transfer without 
reverse tunneling takes place near the negative lead. For 
the system we have studied, the amplitude of this ex- 
cited cycle is about 0.3 A larger than that of the simple 
cycle, and the current level is increased by approximately 
20% to 0.23 /iA. The possibility of two stable shuttling 
amplitudes has recently been discussed by Usmani and 
co-workers 44 



IV. DISCUSSION 

Both the energy spectra and the (Kohn-Sham) or- 
bitals of small molecules near metal surfaces exhibit a 
rich structure and vary substantially as a function of the 
molecule-metal surface separation. The transition rates 
are largely exponential functions of the tunneling dis- 
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FIG. 8: The average current, nn/t, where njj is the number 
of electrons from the negative lead. Inset: The shuttle period 
as a function of time. The shuttle quickly reaches a stable 
motion. 
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tions are forbidden by energy considerations. This results 
in an asymmetry in possible charge states and in asym- 
metric shuttling where energy is absorbed only during 
one half cycle of the periodic motion. 

The forces in the system are highly sensitive to the de- 
scription of the electrode-molecule interaction and to the 
electronic structure of the shuttling object. In the sta- 
ble shuttling regime the island velocity is large enough for 
the island to bounce between the repulsion walls, and the 
details of the forces near the middle of the gap are less im- 
portant than the balance between dissipation and short 
range surface forces. The short range forces are hard to 
describe quantitatively due to the increased importance 
of many-body effects, deformation of molecular orbitals 
near surfaces and the details of the dynamic charge trans- 
fer processes. However, for large enough speed on impact 
the mobile molecule may bounce off the surface and es- 
tablish stable periodic motion. 

The shuttle excitations depicted in Fig. [§l are an ex- 
ample of effects that arise due to the details of the en- 
ergy spectrum of a small system. For a more compli- 
cated spectrum and larger bias voltage more phenomena 
of the same type can be expected; for a slightly different 
model of shuttle-surface interactions we have even ob- 
served that the regular shuttling motion may pass into 
a more chaotic behavior. Therefore, it is likely that a 
microscopic picture of both forces and energy levels is 
paramount for both quantitative and qualitative predic- 
tions of molcculc-sizcd shuttles. 



FIG. 9: The average current for a path with system exci- 
tation due to back-tunneling close to the negative lead (see 
inset). The increase of amplitude caused by the larger Q = 1 
forces enhances the possibility for the tunneling triplet to take 
place also in the next period. The transition is sharp with an 
amplitude increase of ~ 0.3 A and a current of ~ 0.23 fj,A. 
The excitation lasts until a period without reverse tunneling 
takes place. 

tance as assumed in phenomenological theories, but the 
allowed transitions are determined by the energy spec- 
trum, and in particular near the surfaces certain transi- 
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